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Abstract — An emerging and challenging area in mathemat- 
ical control theory called Ensemble Control encompasses a 
class of problems that involves the guidance of an uncountably 
infinite collection of structurally identical dynamical systems, 
which are indexed by a parameter set, by appljing the same 
open-loop control. The subject originates from the study of 
complex spin dynamics in Nuclear Magnetic Resonance (NMR) 
spectroscopy and imaging (MRI). A fundamental question con- 
cerns ensemble controllability, which determines the existence 
of controls that transfer the system between desired initial and 
target states. For ensembles of finite-dimensional time-varying 
linear systems, the necessary and sufficient controllability con- 
ditions and analytical optimal control laws have been shown to 
depend on the singular system of the operator characterizing 
the system dynamics. Because analytical solutions are available 
only in the simplest cases, there is a need to develop numerical 
methods for synthesizing these controls. We introduce a direct, 
accurate, and computationally efficient algorithm based on 
the singular value decomposition (SVD) that approximates 
ensemble controls of minimum norm for such systems. This 
method enables the application of ensemble control to engi- 
neering problems involving complex, time-varying, and high- 
dimensional linear dynamic systems. 

I. Introduction 

The implementation of all scientific and engineering appli- 
cations is complicated by uncertainty or variation in system 
model parameters, for which known control techniques are 
unable to successfully compensate. This issue is especially 
challenging when the control task must be accomplished 
without feedback, whether the control function must transfer 
a single control system between states of interest without 
sensitivity to an uncertain parameter set, or steer a possibly 
uncountable collection of structurally identical systems with 
variation in common parameters between states that may 
depend on the parameters. Investigation of the latter category 
is motivated by factors that arise in practical apphcations of 
quantum control theory to the fields of Nuclear Magnetic 
Resonance (NMR) spectroscopy and imaging (MRI), and has 
given rise to a new area of mathematical control theory called 
ensemble control [1]. Rapidly progressing technologies based 
on quantum theory require the manipulation of very large 
ensembles of quantum systems on the order of Avogodro's 
number (6 x 10^^), whose states cannot be measured during 
the transfer, and whose dynamics are subject to dispersion 
in parameters such as frequency. The performance of the 
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necessary controls must be insensitive to parameter variation 
across the ensemble, as well as to inhomogeneity in the 
applied radiofrequency (RF) control field [2], [3J. A long 
standing problem of significance to NMR requires the design 
of RF excitations that steer a given quantum ensemble 
between initial and target states, and whose performance is 
insensitive to parameter variation [4], [5], [6]. An acceptable 
control function must concurrently drive a collection of sys- 
tems, with identical dynamics but parameter values unknown 
up to a given range, between desired initial and target states. 

Although first motivated by the necessity to control large 
collections of similar systems, the mathematical devices 
produced by investigating ensemble control can also be used 
to approach any open-loop control application in which the 
system response must be immune to uncertainty in model 
parameters. For instance, harmonic oscillators are widely 
used to approximate periodic phenomena in a variety of 
scientific and engineering applications where the frequency 
of oscillation may not be known exactly, but rather is 
confined to a given range [7], [8]. Harmonic oscillators 
often appear in quantum-electrodynamics, and steering such 
quantum systems using electromagnetic fields is a subject of 
widespread interest [9]. 

The theoretical investigation of ensemble control begins 
with the notion of ensemble controllability, which deter- 
mines the existence of controls that achieve various types 
of state transfers for a system of interest. It has been 
shown that a bilinear system evolving on SO (3) called the 
Bloch equations, which models the evolution over time of 
a sample of nuclear spins, is ensemble controllable [3], 
[1]. In addition, the necessary and sufficient conditions for 
ensemble controUabihty of finite-dimensional time-varying 
hnear systems for transfers between states in Hilbert space 
have been recently derived [10]. These conditions depend 
on the singular system of the operator characterizing the 
system dynamics, which can in turn be used to represent 
the minimum norm control that accomplishes the transfer as 
an infinite sum of weighted eigenfunctions. This method was 
used to synthesize optimal ensemble controls for a harmonic 
oscillator system, for which the resulting eigenfunctions are 
the well-known family of prolate spheroidal wave func- 
tions [11]. This special structure facilitates synthesis of the 
controls in this special case, as well as the computation 
of optimal controls with bounded amphtude by solving a 
constrained convex optimization problem [11]. 

The controllability conditions for general, possibly non- 
hnear, ensemble control problems are presently unknown, 
and generalized analytical control design methods remain 



a challenging problem, although analytical solutions exist 
for a few specific systems. When the singular system of 
the relevant integral operator is unavailable in analytical 
form, or in the case of a nonlinear ensemble system, an 
alternative is to use a pseudospectral numerical method that 
translates an optimal control problem in function space into 
a discrete nonlinear programming problem [12], [13], [14], 
[15], [16]. This method has been effective for solving a 
variety of ensemble control problems, but may be difficult 
to implement for large-scale systems with variation in many 
parameters because of the computational complexity required 
for approximating control functions with sufficient accuracy 
[17]. Therefore, a need exists for direct and computationally 
efficient numerical methods for synthesizing ensemble con- 
trols that can accomplish various state transfers for a variety 
of systems. 

In this paper, we introduce an accurate, stable, and com- 
putationally efficient numerical method based on the sin- 
gular value decomposition for constructing minimum norm 
ensemble controls for finite-dimensional time-varying linear 
systems. In Section |ll] we briefly review the theoretical work 
that forms the foundation of our approach. In Section III 
we describe our method for numerically approximating the 
singular system of the Fredholm integral operator of the first 
kind that characterizes the dynamics of a linear ensemble 
system, as well as the synthesis of the unique minimum norm 
control that accomplishes a desired transfer in function space. 
In Section |IV] we revisit the control of harmonic oscillators 
to demonstrate the effectiveness of our new method for 
accomplishing complex state transfers, and also examine 
ensemble control of linear systems with higher dimension 
and with variation in several parameters. Finally in Section 
[V] we conclude by discussing the advantages of the method 
presented here, as well as our future work on developing 
fast, iterative methods for ensemble control of nonlinear 
systems. Such techniques will accelerate the rapidly expand- 
ing scope of ensemble control theory and spectral methods 
by contributing powerful new tools for solving cutting-edge 
problems in diverse fields from neuroscience to quantum 
physics. 

II. Ensemble control of linear systems 

The aim of ensemble control is to simultaneously manipu- 
late a continuum of dynamical systems, which are governed 
by internal and external dynamics that depend on a parameter 
varying over a compact indexing set, by applying the same 
open-loop control input to each. In this section, we review 
the basic definitions and the fundamental theoretical results 
that enable ensemble control synthesis for finite-dimensional 
time- varying linear systems. 

Consider a parameterized family of dynamical systems 
indexed by a parameter f3 varying over a compact set K, 
given by 

X{t, /3) = A{t, (i)X{t, P) + B{t, P)u{t), (1) 
XgMcM", ;3eXcM, tteJ/cM", 



where A{t,l3) e M"^" and e M"^™ have elements 

that are real Loo and L2 functions, respectively, defined on a 
compact set D = [0, T] x K, and are denoted A e L'^'\D) 
and B G L2^™{D). The ensemble controllability conditions 
for the system ([TJ depend on the existence of an open-loop 
control u : [0,T] — ?> U that can steer the instantaneous 
state of the ensemble X{t, ■) : K ^ M between any 
points of interest in the Hilbert space of functions on K. 
Let T-Lt = L2''[0,T] denote the set of m-tuples, whose 
elements are complex vector-valued square-integrable mea- 
surable functions defined on < t < T, with an inner 
product defined by 

{9,h}T= r g\t)h{t)dt, (2) 



where f denotes the conjugate transpose. Similarly, let T-Lk 
L2 (^) be equipped with an inner product 



K 



pHP)q{P)dtl{f3), 



(3) 



K 



where /i is the Lebesgue measure. With well-defined addition 
and scalar multiplication, Ht and Hk are separable Hilbert 
spaces, where 1 1 • | |t and \ \-\\k denote their respective induced 
norms. 

Definition 1: (Ensemble controllability [10]) We say that 
the family ([T]) is ensemble controllable on the function space 
Hk if for all e > 0, and all Xo, Xp G Hk, there exists T > 
and an open loop piecewise-continuous control u £ Ht, 
such that starting from X(0,f3) = Xo{f3), the final state 
X{T,P) ^XtW) satisfies \\Xt ~ Xf\\k < £■ 
In other words, the system ([T]i is ensemble controllable if it 
is possible to guide it from X^ to Xp in the space Hk, 
where the acceptable range of T e (0, 00) may depend 
on e, K, and U. Necessary and sufficient conditions have 
been determined for the ensemble controllability of finite- 
dimensional time-varying linear systems, and are based on 
the Fredholm integral operator that characterizes the system 
dynamics [11]. Given the initial state X{0,(3) = Xo{/3) of 
the system ([T]i, the variation of parameters formula gives rise 
to the solution 

X{T,f3) = ^T,0,p)Xo{l3) 

+ f ^T,a,P)B{a,P)u{a)da, (4) 
Jo 

where $(r, 0, /3) is the transition matrix for the system 
X{t,/3) = A{t,l3)X{t,l3). Our goal is for the terminal 
state to equal the target state in the function space Hk, 
so setting X{T,l3) = XpiP), pre-multiplying by $(0,T,;9) 
and rearranging results in the integral operator equation 

[Lu){fi)^ ( ^{Q,a,P)B{a,P)u{a)dcj^m. (5) 



where = $(0, T, /3)Xf(/3) - Xq{I3). The theory of 

ensemble controllability and the derivation of minimum norm 
controls can be reduced to the solvability of the above inte- 
gral equation. A spectral decomposition, called the singular 
system, of the operator L is used to produce an infinite 



eigenfunction series expansion for the u E Ht of minimum 
norm that satisfies (j5]l with sufficient accuracy. 

Definition 2: Singular System [18]: Let Y and Z be 
Hilbert spaces and L : Y Z he a compact operator 
If (cr,^, t'n) is an eigensystem of LL* and ((T^,/i„) is an 
eigensystem of L*L, namely, LL*Vn = cfn^n, € Z, and 
L*L/i„ = cr^/^n, /!„ e Y , where cr„ > (n > 1), and the 
two systems are related by the equations i/i„ — (7„f„ and 
L*Vn = cr„/i„, we say that (ct„, /i„, Vn) is a singular system 
of L. 

Suppose that ((t„, /x„, is a singular system of the operator 
L as defined in (|5]l, which is compact [10]. The necessary 
and sufficient conditions for ensemble controllability of the 
system ([T]i have been proven to be 



< OO, 



ee7^(L), 



(6) 
(7) 



where 'R{L) denotes the closure of the range space of L. In 
addition, it was shown that the control law 



(8) 



satisfies {u.u)t < {uq,uq)t for all uq € U and mq ^ u, 
where U — {v \ Lv ~ ^}. For further details we refer the 
reader to the original work on this subject [10], [11]. 

Because singular systems and hence optimal ensemble 
controls cannot be derived analytically except for in the 
simplest cases, an accurate and direct numerical method for 
approximating the former is a prerequisite for applying this 
new theory. Given an appropriate numerical approximation 
to the singular system (cTniMm^n) for the operator L of 
an ensemble controllable system, the series ([8]l can be 
truncated to synthesize an approximation to u that results in 
\\Xt — Xf\\k < £ as desired. In addition, a numerical test 
of the criteria (j6]l for ensemble controllability is a natural 
extension of such a framework, which we present in the 
following section. 

III. Numerical approximation of ensemble 

CONTROLS 

The singular system as defined above is the infinite- 
dimensional analogue of the well-known singular value de- 
composition (SVD) for matrices [19]. A natural approach is 
therefore to approximate the action of the compact operator 
L : Ht in equation (j5]l on a function g E Ht by a 

matrix acting on an appropriate vector of sampled values of 
g. Then the SVD can be used to approximate the singular 
system of the operator, and thereby also the minimum norm 
ensemble control u. Let {Pj} be a finite collection of points 
that are distributed uniformly throughout the space K and 
indexed by j = 0, 1, 2, . . . , P, and let {t^} be a collection 
of points that linearly interpolate the time domain [0, T] for 
k — 0,1, . . . , N, including endpoints, with tk — tk-i — 5. 



Using this grid of nodes, we make the approximation 

[Lgm^ ( ^{Q,t,p)B{t,p)g{t)dt 



Jo 

= J2 H0,t,P)B{t,l3)g{t)dt 

N 

^Y.5<^{QM,P)B{tu,p)g{tk). (9) 

The action of the operator i on a function g E Ht can 
be approximated by the action of a block matrix W E 
jgnPxmJV^ with n X m blocks Wjk = <5$(0, tk, l3j)B{tk, l3j), 
on a vector g E M™^, with N blocks gk = g{tk) of 
dimension m x 1. If the SVD of this matrix is = UT,V\ 
and Uj and Vj are columns of U and V, respectively, 
corresponding to the singular value Sj, then WW^Uj = s'jUj 
and W^Wvk — sjvk- Therefore the SVD {sj,Vj,Uj) of the 
matrix W approximates the singular system {aj , Hj , Vj ) of 
the operator L, where Vj and Uj are discretizations of /ij 
and Vj, respectively. Now suppose that ^ E M"^ is given by 
— i{tk) for a function E Hk- Then the minimum norm 
solution g* that satisfies Wg — ^ is given by g* = z 
where WW^ z = ^ [20], so applying basic properties of the 
SVD results in 



E 



(10) 



The components of the synthesized minimum norm control 
u* = {u\, . . . , u^)^ are therefore given by 



^ Sk+,n{j-l) 



'^k+m{j-l)- 



(11) 



Note that the time and parameter discretizations N and P 
must be chosen such that nP < mN, so that the pair 
{W, ^ ) represents an underdetermined system and therefore a 
minimum norm and not a least squares problem. The number 
q of eigenfunctions used in the approximation is limited by 
q<P. 

The use of a Riemann sum quadrature formula to approx- 
imate the action of a Fredholm integral operator of the first 
type by a matrix, so that the SVD can be used to approximate 
the singular system of the operator, has previously been 
examined as part of a least-squares type method [21]. An 
analysis of numerical methods for approximating solutions to 
Fredholm integrals of the first type, as well as an examination 
of accuracy and conditioning issues, has also been performed 
[22]. The latter work includes a discussion of the Picard 
criterion, which asserts that there exists a square integrable 
solution to the integral equation (jSj) only if (i) holds in 
(j6]l. A check of the Picard criterion might provide a test 
for ensemble controllability, but unfortunately it is very 
problematic to test this asymptotic property numerically. 

The most important computational issue is preventing the 
aggregation of numerical errors. These first arise from com- 
putation of the flow which must be done using numerical 
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Fig. 1. Simulation of system ensemble (12) for TV = 20000, T = 1, 
P = 20, and uj £ [—10, 10]. Tlie initial and target states are Xo{uj) = 
(l,0)t and Xp{lj) = (0,0)^. The matrix W is computed in about 5 
seconds, and the SVD is computed in under 1 second, (a) The optimal 
control law {u{t),v{t)) for t £ [0,1] (left), and the final states for all 
systems ui £ [—10, 10]. (b) The singular values {sj} of on a logj^Q 
scale, with the si/sj < 10* cutoff indicated. Here q = 9 singular vectors 
ai'e used to synthesize the control. 



integration when the system is time-varying. A practical 
relative error tolerance for solving ODE systems is 0(10^^). 
Another source of numerical error arises from computation of 
the SVD. We have found that in order to prevent these errors 
from dominating the synthesized control, it is appropriate to 



choose q in (lOi such that the corresponding first and last 
singular values used satisfy si/smq < 10'*. 

IV. Examples and Discussion 

In order to illustrate the performance of our method, we 
will present several examples in this section. The com- 
putations are performed using IMATLAB on a PC with a 
3.33GHz processor. We first revisit the optimal control of an 
ensemble of harmonic oscillators, which has been previously 
examined in detail [11]. This system was proven to be 
ensemble controllable, and the eigenfunctions of the operator 
that characterizes the system dynamics are related to the 
family of prolate spheroidal wave functions. The dynamics 
are given by 



(12) 




-4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 -2.i 
(a) '°9io(8) 

Number of significant singular values q vs. T 




Fig. 2. The simulation in Figure [T] is repeated for different values of time 
horizon T and time step 5, where N = T/5, and P = 40 is used in each 
instance, (a) The norm of the error in the final state is plotted as a function 
of 1/(5 = N/T. The slope of the lines is very close to 1 , so that the error 
is proportional to 5. The lines coiTespond to T = 0.1, 0.5, 2, 1, and 5, 
from top to bottom, hence a longer time horizon does not necessarily result 
in improvement, (b) The number of significant singular values q is plotted 
as a function of T. 



where cj e K = [aji,cj2] C M, the instantaneous state is 
= (a;(-, oj), y(-, o;))^ € T-Lk, and the control vector 
is [/ = (ui,M2)^ € T-Lt- We apply the method described in 
Section [III] to solve an optimal ensemble control problem for 
the system ( 12 1, and the results are shown in Figure [T] A 4th 
order Runge-Kutta scheme is used to integrate the system 
using a relative error tolerance of 0(10^^) for adaptive 
time stepping, and the synthesized control ii* is linearly 
interpolated at the time points requested by the routine. The 
control is similar in shape but not equivalent to previous 
results [11], and its performance is similar. In addition, the 
results of further numerical experiments that determine the 
sensitivity of the error \\Xt — XpWj^ on the choice of time 
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horizon T and time discretization 5 = T/N are shown in 
Figure [2] The error in the final state of the ensemble is 
proportional to the time step 6, which provides a means to 
calibrate the number N of time discretization points. 

This method can also be used to solve more challenging 
problems, for example where the initial and target states Xq 
and Xp are functions of u) in the state space. We applied 
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Fig. 3. Simulation of system ensemble (12) for Af = 20000, T = 40, 
P = 89, and ui £ [—10, 10], where the initial target states Xq and Xp are 
aiTangements of the P -I- 1 oscillators in star' and leaf shaped images in the 
plane, respectively. The average error in the final states of the oscillators 
is 3.03 X 10~*, and the maximum error is 0.0167. (a) Initial and actual 
final states are plotted, (b) The control that accomplishes the transfer, (c) 
The spectrum of the SVD is shown. Observe that the spectrum of the SVD 
differs in form from that which results from the simulation in Figure ^ 
Because the conditioning criterion s\/sq < 10^ is satisfied, all of the 
singular vectors are used to synthesize the control. 



our approach to steer an ensemble of harmonic oscillators 



(12 1 from an initial state arranged in the shape of a star 
to a target state in the shape of a maple leaf in the plane. 
The results of this simulation are shown in Figure [3] and 
we encourage the reader to view a video of the transfer that 
is available online [23]. This example is in fact related to 
a complex problem of importance to the field of NMR in 
which the initial and target states also depend on system 
parameters. In certain experiments specific sub-collections of 
quantum systems must be excited based on parameter values 
or the physical position in the sample under study by using 
so-called selective pulses [5], [24]. Controls that can create 
arbitrary patterns in the terminal state of the ensemble as a 
function of system parameters are therefore desired. 

System ensembles governed by dynamics of higher di- 
mension with variation in several parameters can also be 
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Fig. 4. Simulation of the ensemble system for N = 10000, 

T = 1, r G [-.01,. 01], c e [-.l,.l], andF = 104. The initial 
and target states are xo Ri (0.83,1.38,-1.06,-0.47)+ and xp = 
(-0.27, 1.10, -0.28, 0.70)t. The en'or between the terminal and target 
states is \\Xx — Xp\\fc ~ 0.038. There are mq = 12 eigenvectors used 
to synthesize the control, (a) The minimum norm control, (b) Trajectories 
in the first 3 coordinates, (c) The spectrum of the SVD decays quickly, as 
in the example in Figure [T] 



controlled using our approach. Consider for example the 
time-varying system 

X(t, r, c) = [Ao + Ai sin(27rt) + A2r]X{t, r, c) 

+ [Bo + Bij^^+B2c]u{t), (13) 

where X{t,r,c) e R'^, u{t) e IR^, r e [ri,r2] and 
c e [ci,C2] are parameters, and Ao,Ai,A2 G M.'^^'^ and 
Bo,Bi,B2 E M.'*'^^ have random entries generated using 
a standard normal distribution. Suppose that the initial and 
target states Xo{r,c) = xq E R'^ and Xi?(r, c) ^ xp E 
also similarly randomly generated. Although ensemble 
controllability properties of this system are not straightfor- 



ward to determine, a control synthesis can nevertheless be 
attempted, with a successful outcome indicating ensemble 
reachability of the target state from the initial state. In Figure 
|4]we provide the results of a simulation where a randomly 
generated system of the form ( [T3] l is driven between two 
randomly generated states. This system is quite sensitive 
to variation in the parameter r, while the algorithm can 
compensate for dispersion in the parameter c well. 

The theoretical treatment in previous work provides a solid 
foundation for examirung ensemble controllability [11], but 
a straightforward test for this property is not yet available. 
While it is possible to test for ensemble controllability in 
certain cases by using Lie algebras [2], the complexity of the 
systems encountered in many applications makes this prob- 
lematic in general. It is important to explore the relationship 
between ensemble controllability, the Picard criterion, and 
the singular values of the integral operator (jSj) and its matrix 
analogue W. Investigation in this direction may lead to an 
implementable numerical test for general ensemble systems. 
In addition, a thorough numerical analysis is required to 
better understand the accuracy and conditioning properties 
of this approach, in order to predict its performance under 
various circumstances, and to determine whether a compu- 
tational test for ensemble controllabihty is indeed feasible. 

V. Conclusions 

We have introduced an accurate, stable, and computation- 
ally efficient numerical method for synthesizing minimum 
norm ensemble controls for finite-dimensional time-varying 
linear systems. By basing our approach on the singular 
value decomposition (SVD) and discarding all but the most 
significant singular values, we have guaranteed accuracy and 
numerical stability of our method, and have leveraged the 
efficiency of widely-used numerical routines. Furthermore, 
because the SVD is a finite algorithm, our method does 
not require any additional optimization steps. We have 
demonstrated its effectiveness for designing controls for a 
variety of system ensembles and state transfers under various 
challenging conditions, including complicated state transfers 
and high-dimensional time-varying dynamics. In addition, 
we have conducted multiple simulations to illustrate the 
sensitivity of the method with respect to the chosen time 
step, and the effect of the time horizon on the condition of 
the problem. 

This work forms the basis for a new set of numerical 
methods for quantum control by providing an approach to 
designing excitations for guidance of quantum harmonic 
oscillators. In our future work, we plan to consider fixed 
point and contractive properties of integral operator equations 
to design fast, iterative methods for ensemble control of 
nonlinear systems. We will focus in particular on dynamic 
equations of bilinear form, which govern many widely stud- 
ied quantum dynamical phenomena. A computationally effi- 
cient numerical scheme to synthesize controls for this type of 
problem will facilitate the improvement of pulse design for 
NMR, and would be of immediate use to experimentalists. 
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